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ABSTRACT 


Important properties of derivative (difference) filters 
using the discrete Pourier transform are investigated. The 
filters are designed using the derivative theorem of Pourier 
analysis. 

Because physical data are generally degraded by noise, 
the derivative filter is modified to diminish the effects of 
the noise, especially the noise amplification which normally 
ocours while differencing. The basis for these 
modifications is the reduction of those Pourier components 
for which the noise most dominates the data. 

The various filters are tested by applying them to find 
differences of two-dimensional data to which various amounts 
of signal dependent noise, as measured by a root mean square 
value, have been added. The modifications, circular and 
square ideal low-pass filters and a cut-off pyramid filter, 
are all found to reduce noise in the derivative without 
significantly degrading the result. And the last also 
reduces Gibbs oscillations for those data sets for which 
these oscillations are present with low-pass filtering. 

The FORTRAN programs which perform the filtering 
(DERIV4.F0R and FILT.FOR) and the program which adds the 
noise (GNOISE.FOR) are given and discussed. 
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Filtering Is the process of multiplying the Pourier 
transform of data witn some function. The mathematical form 
of the multiplying function depends on the desired outcome. 
One type of filter that is used to locate peaks , to 
determine the position of boundaries and locate edges, (and 
to determine the derivative) is the derivative filter. The 
drawback to using this filter is the sensitivity of the 
derivative operation to noise, especially high frequency 
noise. This problem can be alleviated by modifying the 
derivative filter. 

In this study a discrete approximation to the 
derivative of the input data is obtained by operating in the 
transform domain using the derivative theorem. For data 
with noise the transform is filtered a second time to reduce 
the effects of noise. Results of these operations with and 
without noise are compared in the function domain. Computer 
programs are used to perform these manipulations of the 
data. 


DERIV4.F0R is the FORTRAN program that outputs the 
transform of the input data (two-dimensional data), the 
transform of the x-derivative, the y-derivative, or the 
second derivative with respect to x and y. DERIV4.F0R uses 
the one-dimensional FFT program listed in Higgins’ article 
(Higgins, 1 976 ) . And DERI V4. FOR accepts real (one input 
file) input data or complex (two input files) input data, 
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while it outputs complex data. The program that adds the 
Gaussian noise to the data is GNOISB.FOR, and it accepts 
only real data. 0N0I3E.P0R is based on a noise program 
written by Bill Bivens (Bivens ,1976) which has been modified 
for two-dimensional data. PILT.FOR is the FORTRAN program 
that performs the filtering, and was written by the author. 
The three filter choices are a circular filter, a square 
filter (called a rect filter), and a flat-topped pyramid 
filter. The filter consists of the multiplication of one of 
three filter functions after the transform has been 
multiplied by the derivative filter. The circular and rect 
filter functions have the value of one inside the boundaries 
of the respective geometrio figures after which the filters 
are named, and have the value zero outside the boundaries. 
The flat-topped filter function is one inside the square and 
slopes linearly (with negative slope) from the edges of the 
square to the end of the data matrix where the value of the 
function is zero. FILT.FOR also allows for the size of the 
filter function (the size of the square or circle) to be 
varied. 

This thesis is separated into three chapters with a 
summary and an appendix. The three chapters are entitled 
FOURIER TRANSFORMS, TAKING the DERIVATIVE, and NOISE and 
FILTERS. The appendix includes the FORTRAN code for the 
programs that were used. 
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CHAPTER 1 

FOURIER TRANSFORMS 


1 1 SPECIAL FUNCTIONS AND THEOREMS 

In one dimension the Fourier transform of a oontlnuoue 
function, f(x), ie t 



K 4 

e 


dx 


This Integral ie a function of s and will be called F(s). 
The independent variables In the function and transform 
domains are x and s, respectively. Here the term frequency 
will be used fco represent the independent variable in the 
transform domain, regardless of the units of the function 
domain independent variable. The inverse transform of F(s) 
is : 


fo - / F(. 
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at?** w 




The sign reversal in the exponential is necessary to ensure 
that two successive transformations result in the original 
function (Bracewell, 1965)* 


3 


4 


Ths Fourier transform of a continuous funotion In two 
dimensions is: 

r . / - /* S (*< ♦ *y> 

F(*,*r)» £ J « J* cfj 

where u Is the independent variable in the transfers domain 
associated with the independent variable x of the funotion 
domain and v of the transform domain is associated with y of 
the function domain. The inverse transform is: 

n*»y) * j J P( a »* r ) c ^ 


First I will define some common functions associated 
with Fourier transforms and give their transforms 
(Bracewell, 1965 ). 

In one dimension some common functions are 


rect(x) 

- 1 

abs(x) 

< 1/2 


« 0 

abs(x) 

> 1/2 


« 1/2 

abs(x) 

* 1/2 

tri(x) 

* 1 - abs(x) 

abs(x) 

< 1 


* 0 

abs(x) 

> 1 


thm Impulse function 

d(x)i 1> / d(x)dx j d(x) b 0, x 0 0 

tho nine function 
sino(x) b 

tho replicating or shah function 
III(x) b J^d(x-n) 

even impulse pair 

Il(x) « l/2d(x+l/2) + 1/2d(x-1/2) 
odd impulse pair 

Ii(x) b l/2d(x+l/2) - l/2d(x-1/2) 


Gaussian 
exp( -ex 4 ) 


COMMON TRANSFORM PAIRS 


f(x) 

rect(x) 

tri(x) 

d(x) 

III(x) 

Il(x) 


P(8) 

sinc(s) 

sinc*(a) 

1 

III(o) 
cos (ft's) 


taints) 

exp(-flb*) 
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Ii(x) 
txp(-^x 4 ) 


In two dimensions the common functions are 


rect(x,y)* 

sinc(x,y)« 

d(x,y)« 

gaussian» 


rect(x)rect(y) 
einc(x)einc(y) 
d(x)d(y) 
exp(-7T(x a + y*)) 


COMMON TRANSPORM PAIRS 
f(x,y) P(u,v) 

rect(x,y) sinc(u,v) 

exp(-7?tx 4 + y*)) exp(-^(u A + v a )) 

C08(^x) Il(u)d(v) 

A few theorems play a basic role when dealing with 
Fourier transforms. Therefore I will state these theorems 
here for those unfamiliar with Fourier transforms, but will 
not give the proofs, which may be found in Bracewell, 1965* 

A few of the basic theorems: 


THE SIMILARITY THEOREM 

If f(x) has the Fourier transform F(s) then f(ax) has 

•I 

the Fourier transform |a|F(s/a). Conceptually this states 
that broadening a function in the function domain causes 
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contraction of the transform and growth of its ordinate in 
the transform domain, and vice versa. 

THE ADDITION THEOREM 

If f(x) and g(x) have the Fourier transforms F(s) and 
Q(s), respectively, then f(x) + g(x) has the Fourier 
transform F(s) + 0(o). 


THE SHIFT THEOREM 
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If f(x) has the Fourier transform F(s) then f(x-a) has 
the Fourier transform exp(-2^ias)F(s) . This simply states 
that shifting a function in the function domain is analogous 
to giving the transform a frequency dependent phase shift. 


THE CONVOLUTION THEOREM 

r 

If h(x) * / f(u)g(x-u)dx (h*f*g, f convolved with g) 
and f(x) has Fourier transform F(s) and g(x) has Fourier 
transform C(s), then H(s) is the transform of h(x) where 
H(s) * F(s)G(s). The convolution theorem is used quite often 
since it is usually easier to perform a multiplication or 
division than it is to perform a convolution or 
deconvolution. 

THE DERIVATIVE THEOREM 
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If f(x) has the Fourier transform F(s) then f*(x),the 
derivative of f(x) hae the Fourier transform 27fisF(s). Here 
is a method for obtaining the derivative of a function that 
is not an approximation technique. And by successively 
applying this theorem higher order derivatives can be 
obtained. 

Another theorem, although not a basic theorem in 
continuous Fourier transform analysis, is important in going 
from continuous transforms to discrete transforms. This 
theorem, known as the sampling theorem, states that a 
function whose transform is zero for s > s 6 (where s c is the 
cutoff frequency in the transform domain) is fully specified 
by samples taken at equal intervals not exceeding l/(2s c ) 
save for any harmonic terms with zeroes at the sampling 
points. Thus it is possible to reconstruct a function from 
its samples if the sampling interval is less than or equal 
to l/(2s e ) . See Fig. 1.1 for graphical detail. 
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'’’ho shah function 1m useful In representing sampled 
data. For wh*>n ono multipliers a function f(x), by the shah 
function, lll(x), one la effectively sampling that 
continuous function at evenly spaced intervals. The values 
of the function, f(x), at integral values of x are preserved 
(by the deltas) whereas information of the function between 
the intervals of the deltas is not kept. Symbolically this 
sampling property of the shah function is represented as: 

Pi 
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And this ia graphically stated In Pig. 1.2. 

Another important property of the shah function that la 
closely associated with the sampling property la that of 
replication. When the shah function is convolved with 
another function the result is the function being replicated 
at unit intervals to infinity in both directions. If the 
function is wider than one unit interval (or wider than the 
replication interval when this is not unity) then there is 
overlapping of the replications. When this occurs in the 
transform domain it is known as aliasing and can be a 
serious problem. The replicating property symbolically 
stated is: 

* P(A = 2 

A* -to 

The replicating property of the shah is shown graphically in 
Pig. 1.3. 

Discrete functions (data) can be viewed as sampled 
continuous functions. That is, if the continuous function 
is represented by f(x) then the discrete (or sampled) 
version can be represented by IIl(x/t)f(x) , where t is the 
sampling interval. And as a consequence of sampling, the 
transform of a discrete function is |t|*III(ts)*F(s) , where 







12 


F(s) 1 b the transform of the oontlnuous function. Thus wtt 
see that sampling in the function domain causes replication 
in the transform domain. This oan also be applied in the 
reverse direction- sampling in the transform domain causes 
replication in the function domain. 

1.2 THE DISCRETE FOURIER TRANSFORM 


The Fourier transform as previously stated operates on 
continuous functions, whereas physical data are normally 
discrete. Therefore one must shift gears and begin thinking 
in terms of discrete functions. One method is to construct 
discrete functions from continuous functions using the 
sampling function, and then to use the continuous Fourier 
transform to obtain a representation for the discrete 
Fourier transform (known as the DFT). 

Thus a discrete function can be represented by 
III(x/t)f(x), where f(x) is a continuous function that 
corresponds to the discrete function. Also since the 
discrete function is taken to have finite extent one can 
represent it as the continuous function times a rect 
function. The rect function is often called a window since 
when it multiplies a continuous function only that portion 
of the function which falls under the rect is left non-zero. 
Thus the discrete function can be represented by the product 
III(x/t)rsct((x-a)/c)f(x) . The parameter t in the shah 


function determines the sampling Interval, the parameter c 
in the reot function determines the width of the window, and 
a determines the position of the oenter of the reot. Prom 
the convolution theorem the transform of a discrete funotlon 
can be viewed an |tc|»exp(-2iyias)ni(ts)*sinc(cs)*P(s) . Thus 
the DPT of a discrete function is not exactly equal to a 
sampling of the Pourier transform of the analogous 
continuous function. Due to the processes of sampling and 
windowing the original function is altered somewhat and thus 
its transform is affected also (Bracewell ,1965) • 

Convolving the continuous Pourier transform with a sine 
broadens it, and the convolution with the shah function 
causes replication. If any portion of a function that is 
replicated extends beyond the cutoff frequency then these 
high frequencies mask as lower frequencies and aliasing 
occurs. Since convolving with the sine broadens the 
transform this contributes further to aliasing. Convolution 
with the sine also causes Gibbs oscillations about any rapid 
change in the transform ( Bracewell, 1 965 ) • 

The following is a derivation of the DPT using the 
representation of discrete functions by continuous functions 
multiplied with the shah function (loup,1978). 
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The discrete Fourier transform of discrete data is a 
sampled function. The transform is sampled to allow 
representation on, and calculation using a digital computer. 
Thus the data in the function domain are replicated. The 
fast Fourier transform technique is Just a quick method for 
obtaining the DPT of a function. Therefore the fast Fourier 
transform (FFT) and the DFT (both are the name transform- 
one is Just a faster approach to calculation for large data 
sets) view the data from the function domain as one period 
of the replication with the period of replication equal to 
the width of the window. 

As a result of this if the period of periodic input 
data is incommensurate with that of the window, then it is 
better to add zeroes to the end of one period of the 
function till the periodicity of the function plus the 
appended zeroes is commensurate with the window. Function 
domain replication, implicit in the use of a sampled 
transform, Joins incomplete periods of the function together 
causing discontinuities to be generated when the data is 
incommensurate with the window. See Fig. 1.4. 
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Pig. 1.4 Incommensurate function without and with seroes 


The discontinuities could introduce Gibbs oscillations, in a 
DPT representation These will be discussed in more detail 
in the chapter on noise and filters. 


The cutoff frequency in the transform domain is equal 
to one-half the Inverse of the sampling Interval in the 
function domain (i.e., »l/(2t)). Thus to Include higher 
frequencies and reduce aliasing the sampling interval is 
made as small as possible. 


1.3 THE PAST FOURIER TRANSFORM 


The FFT algorithm reduces the number of operations 
performed in the calculation of the DPT of a sequence. The 
algorithm was rediscovered by Cooley and Tukey in 1964 
(Cochran et al). The significance of thi9 algorithm is that 
it reduces the time required to calculate the DPT of a 
sequence. It takes N multiplications to compute the DPT in 
the straightforward method, whereas the number of 
multiplications performed using the FFT algorithm is 
approximately 2Nlog*N. As N gets large the savings in 
computation time becomes great. For one-dimensional data 
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M*1024 la commons for thla exaaple tha savings amount to a 
faotor of ona hundred reduction In tha number of oparatlona 
required (Higgins ,1976; Coohran at al,1967). 

Basically tha PPT algorithm can ba understood by taking 
an N point tranafora and apllttlng it Into two W/2 point 
transforms. Than thaaa two tranaforaa ara aaoh apllt In 
half, and thla process rapaato itself until thara ara N ona 
point tranaforaa. Thara ara othar PPT algoritaa that work 
on sequences of N points (N not prlma) for N not two ralaad 
to an IntagraX power , but baoausa tha raduotlon to N ona 
point tranaforaa oan not ba oomplatad thay ara not as 
efficient* Tha PPT algorithm also usas tha pariodioity of 
tha axponantial function to allmlnata redundant operations. 


To understand this process let A(r) ba tha value of tha 
tranafora of X(k) at tha frequency r*rns, where 
r *0,1 ,2, . . .N-1 . N is tha number of points in tha saquanoa 
X(k) . Than from the DPT 


A(r) * jLx(k)exp(-21Urk/N) 

k.« 


Then the data set Is split into even and odd saquenoas, Y(k) 
and Z(k). 


k*0,1 ,2, . . . ,(N/2)-1 . 
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Y(k)« X(2k) 


It 


Z(k)«X(2k+1 ) k-0,1 ,2, (W/2)-1 . 


And let 


A(r)«^CY(k)#xp(-mpk/!l) + Z(k)*xp(-2iri(2k+1)p/M)] 


■ ^ Y(k)«xp(-4^1rk/N) + exp(-2^ir/N) J> j Z(k)exp(-4fl , 'irk/N) 
Whtp# Trn 0,1 ,2, . . . ,H-1 


4*«‘ 


NM-t 

if K(p)« 2)Y(k)exp(-4#lrk/lf) 

k*t 


a&;/ 

L(p)« 2;z(k)exp("-4Plrk/l) 

let 


M(r)» exp(-2M.r/W) 


then 


A(r)« K(r) + M(r)L(r) r*0,1 ,2, . . . ,(N/2)-1 . 

Since K(r) and L(r) are periodic in the half interval 
( 0<«r<N/2 ), A(r) can be generated for the second half 
using the values of K(r) and L(r) for the first half. 

A( r+(N/2) )■ K(r) - M(r)L(r) 


The minus sign comes from exp(-2rir/N) as 0<«r<N/2 being 


It 


opposite in sign to exp(-2ftlr/lf) ss N/2<ar<M, (i.e., 
exp(-2*ir/W) ■ -exp(-2rl(r+(N/2 ))/»)). Thus the M point 
trsnsfors has gone tvo lf/2 point transforms. 

The PPT algorithm used (Higgins, 1976) assuaes that the 
first data point is the value of the funotion at the origin , 
and that any values asoooiated with negative x (absolssa) 
are plaeed beyond the funotion value of the last positive x* 
See Pig* 1.5* The data oust either input to the PPT 
prograa in this fora or a portion of the prograa must be 
devoted to rearranging the data into the foraat required by 
the PPT subroutine. The prograa DBRIV.POR which perforae a 
one-dimensional PPT of a sequence does the latter. The 
output of DERIV.FOR is In the sane foraat as the Input. 
Though instead of rearranging the transforn the input data 
is aultlplled by a phase factor (the input data can be 
viewed as having been shifted to the center of the aatrlx). 
"he result is the same as if the rearrangement had been 
performed ( And raws, 1 970 ) . 

When the number of data points Is a power of two one 
runs into diffloulty in representing even functions that 
have their origin sampled. Since the funotion is 
represented by an even number of points the window oan not 
be symmetric about the origin. Thus it is best to have the 
funotion go to sero at both ends within the width of the 
window. If this is not possible then the asymmetry and the 
replication in the function domain whioh results from having 
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a sampled transform aunt be carefully considered. 

In two dimensions the DPT is defined as: 

f ' *•/ 

2x(3,k)exp(-2^i((aj/M) + (nk/N))) 

£7 |*9 

which oan be rewritten as: 

F(m,n)« Id X( j,k)exp(-2*'imj/M)]exp(-2?'ink/N/ 

The term in the brackets is the transform of row (or column) 
k. And the outer sum transforms the columns of the above 
result. Thus a two-dimensional DPT can be obtained by first 
performing a one-dimensional DPT on each row (this amounts 
to M one-dimensional DFT operations) and then executing a 
one dimensional DPT on each column of the matrix of 
transformed rows (this amounts to N one-dimensional DFT 
operations). This is equivalent to the result obtained when 
columns and rows are interchanged in the above procedure. 

The two-dimensional discrete Fourier transform is also 
sampled. The input data are now a matrix of values, where 
the first index corresponds to the y coordinate values and 
the second to the x coordinate values. 

Since the one-dimensional FFT algorithm used in the 
two-dimensional PPT is of the type discussed previously, the 
two-dimensional PPT routine also expects data in a different 
format than what might be expected. In addition, the 
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transform is not arranged as expected. Instead of the 
origin being located at or near the oenter of the transform 
matrix (at point N/2 + 1, N/2 + 1) for an even matrix (where 
N is the number of rows or columns in the matrix) the PPT 
results in the origin being located at the top left corner 
of the matrix with all the low frequencies in the corners 
and the high frequencies in the center. The same idea 
applies to data in the function domain. 8 gs Fig. 1.6. 

The two- dimensional FFT program used, DERIV4.F0R, 
expects data with the origin located at the center and 
rearranges it into the format expected by the 
two-dimensional PPT subroutine. The two-dimensional 
transform is also rearranged such that the origin is at the 
center by the phase multiplication method mentioned 
previously in the section on the one-dimensional PPT. 

Since the number of points is normally even, and in our 
case N where N is a power of two, there is a missing bottom 
row and a rightmost column if data with even symmetry and 
with the origin sampled are used. Since the transform is 
sampled the input data are viewed as one period of the 
replication. The one-dimensional discussion generalizes to 
two-dimensions . 

As in the one-dimensional case sampling in the function 
domain causes replication in the transform domain. Though 
now the replication is in two dimensions with the top row of 
one period adjacent to the bottom row of another period and 
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CHAPTER 2 

TAKING THE DERIVATIVE 

Prom the derivative theorem the transform of the 
derivative of a function ie just 2^is times the transform of 
the function to be differentiated. For higher order 
derivatives the relationship is! c^(f w (x))» (2irln) n F(s) , 
where c/(f(x))» F(s). 

Thus the derivative theorem provides a method for 
obtaining the derivative of continuous and discrete 
functions without approximation other than any approximation 
already made in treating a continuous function discretely. 
This method, like any derivative technique, is sensitive to 
noise, especially since there is no smoothing due to the 
approximations of normally used numerical techniques. But 
if one is already using the FFT and the data are relatively 
noise-free then this method is definitely a viable 
alternative to derivative approximation methods. Also, if 
the data are noisy very effective filters may be used as 
part of the derivative process. 

Since I have concentrated this investigation on 
two-dimensional functions most of the programs operate on 
two-dimensional data with most of the results for two 
dimensional data. When it comes to explaining the theory I 
will use one-dimensional functions for simplicity, reverting 
to two-dimensional functions only to illustrate an important 
point or a veiled implication. 
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The FORTRAN program DERI V4. FOR takes as input 
two-dimensional data in the form of a two-dimensional square 
matrix, with the maximum sise of the matrix being 64 X 64* 
Most of the time though, the sise of the data matrix was 32 
X 32, a compromise between the number of points desired and 
the length of time required to run the program. Using 32 X 
32 matrices DERI V4. FOR ran in approximately one third of the 
time compared to when 64 X 64 matrices were used. 

DERIV4.F0R gives the user four choices of how the 
transform will be manipulated. Once the transform is 
obtained it is multiplied by 2friu, or by 2?Mv, or by -4^ - uv, 
or it is untouched. Thus DERIV4.F0R can give the transform 
of the x-derivative of the input data, the transform of the 
y-derivative of the input data, the transform of the second 
derivative with respect to x and y of the input data, or the 
transform of the data. To obtain the transform of the 
second derivative with respect to x or y DERIV4.F0R is just 
run twice, with the output from the first run being the 
input for the second run. 

Also DERI V4. FOR performs either the rainus-i or the 
plus-i transform. The minus-i transform has a negative i in 
the argument of the exponential whereas the plus-i transform 
has a positive i in the argument of the exponential. 

In multiplying by 2friu or 2fTiv, etc., the transform is 
"centered” (as centered as can be using an even number of 
rows and columns) about the origin. Thus in some cases the 
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sign of the u or v ie negative. That is, although any 
replication could be used, we use the one oentered about the 
origin. 

The sampled functions that were used to check 
DERIV4.F0R were a two-dimensional gauss ian and a cosine 
wave. The gaussian used was exp[((17-I) + ( 1 7 ) ^ )/(4*0)]. 
Figure 1 .7 is the gaussian and Fig. 1 .8 is the x-derivative 
of the gaussian. First DERIV4.F0R was run to obtain the 
minus-i transform of the x-derivative of the gaussian. This 
was then plus-i transformed to obtain the x-derivative of 
the gaussian. The same sequence of events were followed to 
obtain the x-derivative of the cosine wave. The relation 
used for the cosine wave was 1 + cos(2fT(J-17)/l6) where the 
addition of one was to produce non-negative data. Figures 
1.9 and 1.10 are respectively plots of the cosine wave and 
its derivative. 
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CHAPTER 3 
NOISE AND PILTERS 


3.1 NOISE 


Since noise in signals is very common, whether it be 
background noise, instrument noise, or another unwanted 
signal, developing the method of taking the derivative using 
the derivative theorem of Pourler transform analysis would 
not be complete without Including noisy data. The type of 
noise chosen was ordinant dependent Gaussian additive noise, 
i.e., noise with a Gaussian probability density function. 
This choice was made since the noise associated with most 
imaging sensors can be modeled as a Gaussian distributed 
random process. 

GN0I8E.F0R is the FORTRAN program that adds noise to 
the input data. The output is noise added onto the data. 
If f(y) is the probability density function then f(y) 
* e . Now f(y) is the probability density 
function of the noise. To determine the amplitude of the 
noise at a particular point the relation between the 
amplitude and probability density function must be 
determined . 
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Since f(y) is evenly distributed, f(y) can be represented by 
a uniformly distributed random number. Thus, if P is a 
random number between zero and one, y ■ jr-J 1 

To describe ordinant dependent noise «*' is not a 
constant but is equal to a scale factor times the ordinant 
of the data point in question ( *SF»A(I , <J) , where A(I,J) 
is the value of the function associated with the point 
( ( J-17) , (17-1) ) ) . The scale factor allows the root mean 
square value (RMS) of the noise to be varied. 

The amplitude of the noise is y. This is added to the 
ordinant (or the value of the function associated with the 
point) by GNOISE.FOR. Additive noise was used since it is 
common and the simplest to deal with mathematically. 

The synthetic data sequences to which noise was added 
were the gaussian and cosine wave used before* The scale 
factors used were 0.00001 and 0.0001 . The following plots 
show the data with noise and then its derivative. For the 
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Gaussian the derivative ie real, thus the imaginary part 
gives an idea of the round-off error* For example, the 
magnitude of the maximum and minimum values of the imaginary 
part of the x-derivatlve of the Gaussian data are 42*749 and 
-42*349 respectively for the eoale factor equal to 0.00001. 
For the maximum and minimum values associated with the other 
x-derlvatives of the Gaussian data and the oosine wave data 
with the associated scale factors of the noise; see Table 
I. 


TABLE I 


d/dx of 

maximum 

minimum 

scale factor 

Re/lm 

Gauss . 

42.3 

- 42.3 


.00001 

Im 

Gauss . 

12492.0 

-1 2655.5 


.00001 

Re 

Gauss * 

124*5 

- 124.5 


.0001 

Im 

Gauss * 

12018.8 

-13681 .5 


.0001 

Re 

oosine 

274.2 

- 274.2 


.00001 

Im 

cosine 

7589.0 

- 7439.5 


.00001 

Re 

cosine 

772.4 

- 772.4 


.0001 

Im 

cosine 

11306.4 

- 9536.4 


.0001 

Re 

The RMS is 
That is 

the root 

mean square 

of the 

noise 

RMS* 

Zi 

< 

r 

AT 

= rlocse 

i 

2 number o 



The SNR ie the signal to noise ratio, which is 
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SNR« PEAK SIGNAL VALUE 

MIS of NOISE 

Each filter operated on data with the nolee scale 
factor equal to 0.0001 and 0.00001 . Also the SNR and RMS 
for the Gaussian and cosine waves with noise soale faotors 
of 0.0001 and 0.00001 are listed In Table II. 

TABLE II 


function 

soale factor 

SNR 

RMS 


Gauss . 

.00001 

814.3 

K\ 

CM 

• 

B-2 

Gauss . 

.0001 

285.6 

.350 

CM 

S 

m 

cosine 

.00001 

188.1 

.106 

B-1 

cosine 

.0001 

66.7 

.300 

E-1 

3.2 FILTERING 





The aim of any filtering is to reduce the unwanted 
effects in data and enhance the wanted ones. We wish to 
minimize the effects of the noise on the derivative of the 
data. We accomplish this by filtering in the transform 
domain, because if any of the characteristics of the noise 
in the transform domain are known then it is simpler to 
design the filter in the transform domain. Most often the 
only characteristics of the noise that are known are those 
which are simply given in the transform (frequency) domain. 
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In the transform domain ths transform of the dsrlvatlvs 
of ths function is 201s P(s), where P(s) ;'s ths transform of 
ths function itsslf. Ths multiplication by s amplifies high 
frequency noise, therefore any filter should decrease this 
effect. One method would be to out off any frequencies 
above a certain value* This is analogous to multiplying the 
transform of the derivative with a two-dimensional rect 
function, a two-dimensional circular function, or some other 
geometric shaped plateau function. The drawback with this 
type of filtering is the Introduction of Oibbs oscillations 
due to the abrupt windowing in the frequency domain by the 
filter. To reduce any Gibbs oscillations a tail can be 
added to the plateau filter, though this also increases the 
effect of high frequency noise on the transform. 

Gibbs oscillations are the oscillations that result 

around rapid changes in the function domain when the 

function is represented by a transform that has been 

truncated ( multiplied by a rect function in the simplest 

case) (Bracewell,1965) . A truncated transform translates in 

the function domain to convolving the function with a sine 

function, if the region of a discontinuity in the time 

domain is to be examined. The discontinuity can be 

approximated with the sgn function( sgn(x)* 1, x>0; ■ -1 , 

x<0). Now sgn(x)*sinc(x) * 2/sinc(t)dt , where 

* * 

2^sinc(t)dt * (2/p')Si(vrx) (Si is the sine integral). This 

function oscillates about >1 for large negative x values. 

As the origin is approached the amplitude of the 
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oscillations increases, passes though stro at x«0, shoot* up 
to a maximum of 1.18 and th*n oaolllato* about 1 aa x 
Increases, with the oscillation* dying out aa x increases. 
The amplitude of the oscillations about -1 and 1 remains the 
sane if the aino function is compressed by a faotor of N and 
strengthened by a faotor of If ( to preserve unit area) and 
only the frequency of the oscillations ia altered. It ie 
Increased. Thus changing N does not change the amount of 
overshoot which is approximately nine per cent of the amount 
of the discontinuity. 

Thus to reduce Gibbs oscillations a linear tail was 
added to the reot function to form a flat-topped pyramid 
function. This flat-topped pyramid filter function, the 
oircular filter function and the reot filter funotion were 
used on noisy gaussian and cosine wave functions. 

Dealing first with the x derivative of the Gaussian 
data, the worst filters (worst in terms of affecting the 
presence of the noise) were the circular filter with the 
radius, R, equal to four and the rect filter with the length 
of a side equal to nine (*2M+1 , where M»4), with the 
circular filter being worse. (Because the derivative of the 
Gaussian should be real, even with noise, the non-aero 
imaginary parts reflect round-off error in the 
calculations) . 
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For the circular filter with R«4, the oscillations in 
the derivative were quite large and were not dying down 
within the period* While the oaoillatione in the derivative 
froa the reot filter with M«4 were still rather large, they 
were dying off soae as x and y varied froa sero. 

With R«8 the osoillations in the reeult using the 
oiroular filter are less and show soae oiroular symmetry. 
The result with the r»ot filter with M«8 has osoillations 
whioh are lees than for M»4, and die off more rapidly for 
values of x and y off axle. Also the oscillations are 
greatest in the x direction. This is the direction in which 
the derivative is taken, and the derivative process tends to 
amplify the effects of any type of noise. 

The pyramid filter adds a linear ramp to the rect 
function that extends from the edge of replication (or edge 
of the square containing all the data- not including the 
extra row or column) to the edge of the reot. The 
derivatives using the pyramid filter show no perceptable 
oscillations. 

For the cosine wave data the rect and circular filters 
will work the best since the transform of a cosine wave is 
the even pair situated on the u axis with their separation 
determined by the pei iod of the cosine. The transform of 
the cosine wave chosen for this work is non-zero only close 
to the origin. Therefore the rect and circular filters can 
outoff much of the transform ( and thus much of the noise ) 


but still have only a snail effect on the even pair 
(actually the even pair convolved with a slno because of the 
windowing in the function domain). 
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This was seen in the results of the circular, rect, and 
pyramid filters. The pyramid filter was the worst since it 
let in more noise. And as expected as R or M decreased in 
the circular, rect or pyramid filters the results showed 
less effects of noise. 


The cosine wave is representative of functions whose 
spectrum is centered closely about the origin and such 
drastic measures (M*4,R*4) would eliminate important 
information In the transform domain for functions which are 
not so concentrated. Thus, since the Gaussian had a 
spectrum which spread out over the entire period 
(two-dimensional) in the transform domain, the pyramid 
filter was better as it allowed more information of the 
transform through. 
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SUMMARY 

A study of derivative filters using the discrete Pourier 
transform has been performed. As has been discussed a 
filter multiplies the transform of the data to be filtered 
with some function. Thus the derivative filter multiplies 
the transform of the data with 2#'is, where 8 is the 
Independent variable in the transform domain. This result 
is the derivative theorem of Fourier transform analysis. 

But because the derivative process is sensitive to 
noise, the filter must be modified to reduce the noise 
effects. This noise filtering can be done before, after, or 
combined with the derivative filtering since the 
multiplication is commutative. In this case the noise 
filtering was done after the derivative filters. 

The input data were a two-dimensional cosine wave and a 
two-dimensional Gaussian wave. The input data were entirely 
two- dimensional, though in the theoretical portions of this 
thesis reference was made to one-dimensional functions for 
simplicity. And since the data were two-dimensional all the 
computer programs were written to operate on two-dimensional 
data . 

The Gaussian ordinant dependent noise was added to the 
input data by the program GNOISE.FOR. This program uses the 
random number generator in FORTRAN to determine the size of 
the noise. 
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DER1V4*F0R was used to take the FFT of the input. The 
filter program used vas FILT.FOR. 

After a set of noisy data was filtered by DERIV4*F0R 
and FILT.FOR the resultant output was transformed back to 
the function domain in order to examine the derivative of 
the noisy data after filtering. 

Of the three filters, the one employing a tail ( the 
function sloped down to zero rather than abruptly cutting 
off) worked the best on data having a transform not 
concentrated about the origin. Gibbs oscillations in the 
function domain are reduced by such filters and the 
trade-off in letting more noise through to allow the 
transform to go gradually to zero is definitely beneficial. 

For data with transform information centered about the 
origin, the filters that cut off abruptly worked better than 
the filter that employed a tail. This is due to being able 
to get rid of the high frequency noise without destroying 
any important transform information. Thus the appropriate 
filter depends on the type of data. Since the information 
in the transform domain is generally not concentrated about 
the origin, filter functions that are to be used for common 
data need tails. 

Future work in this area might be to use different 
tails such as a Gaussian tail for the circular filter or a 
cosine wave tail for the rect filter. Also the filter could 
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Perspective Plots 


The following are 

plots of the 

final 

results. S? 

is 

the 

scale factor used 

to determine 

the 

amplitude 

of 

the 

noise. Cir, Pyr, 

and Rect 

are 

abbreviations 

for 

the 


circular, flat-topped, and square filter functions, 
respectively. 

1 .7 Gaussian input data 

1 .8 x-derivative of Gaussian data 

1 .9 Cosine wave input data 


1 .10 

x-derivative 

of 

Cosine 

data 



1 .11 

x-derivative 

of 

Gauss . 

SF-1 . OE-4 , 

no filter 

1.12 

x-derivative 

of 

Gauss . 

SF*1 . OE-3 , 

no filter 

1.13 

x-derivative 

of 

Gauss . 

SFsl .OE-4, 

Pyr, 

M*4 

1.14 

x-derivative 

of 

Gauss . 

SF*1 .OE-4, 

Pyr, 

M*4 

1 .15 

x-derivative 

of 

Gauss . 

SF*1 .OE-3, 

Pyr, 

M=4 

1.16 

x-derivative 

of 

Gauss . 

SF*1 .OE-4, 

Pyr, 

M*8 

1.17 

x-derivative 

of 

Gauss . 

SF*1 .OE-3, 

Pyr, 

M-8 

1 .18 

x-derivative 

of 

Gauss . 

SF=1 .OE-4, 

Cir, 

R*4 

1.19 

x-derivative 

of 

Gauss . 

SF=1 .OE-3, 

Cir, 

R=4 

1 .20 

x-derivative 

of 

Gauss . 

SF*1 .OE-4, 

Cir, 

R=8 

1 .21 

x-derivative 

of 

Gauss. 

8P«1 .OE-3, 

Cir, 

Ra8 

1 .22 

x-derivative 

of 

Gauss . 

SF«1 .OE-4, 

Rect 

M*4 

1 .23 

x-derivative 

of 

Gauss . 

SF=1 .OE-3, 

Rect 

M=4 

1 .24 

x-derivative 

of 

Gauss . 

SF=1 .OE-4, 

Rect 

M=8 

1 .25 

x-derivative 

of 

Gauss . 

SF=1 .OE-3, 

Rect 

M=8 
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.1 .26 

x-dorlvntive 

of 

oosine 

3P« 

1 .OB-4, 

Rect 

M»8 

1.27 

x-derivative 

of 

ooeine 

SF« 

1 .OB-3, 

Rect 

M«8 

1.28 

x-derivative 

of 

ooe ine 

SF-1 .OB-4, 

Cir, 

R*8 

1 .29 

x-derivative 

of 

cosine 

SP« 

1 .OB-3, 

Cir, 

R*8 

1 .30 

x-derivative 

of 

cosine 

SP- 

1 .OB-4, 

Cir, 

R*4 

1.31 

x-derivative 

of 

cosine 

SP« 

1 .OB-3, 

Cir, 

R*4 

1 .32 

x-derivative 

of 

cosine 

SF« 

1 .OB-4, 

pyr. 

M«8 

1.33 

x-darivative 

of 

cosine 

SP« 

1 .OE-3, 

Pyr, 

M*8 

1.34 

x-derivative 

of 

cosine 

SF* 

1 .OB-4, 

Pyr, 

M«4 

1.35 

x-derivative 

of 

cosine 

SP« 

1 .OE-3, 

Pyr, 

M*4 

1 .36 

x-derivative 

of 

cosine 

SP* 

1 .OB-4, 

no filter 

1 .37 

x-derivative 

of 

cosine 

SF» 

1 .OE-3, 

no filter 
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Appendix 


Following is the code for the programs discussed in 
this thesis. The order of the execution of these programs 
is s step one- input the data to GNOISE.FOR and choose the 
scale factor of the noise; step two- input the noisy data 
to DERI V4. FOR and choose the minus-i transform and the x, 
second xy, or y derivative; step three- filter the output 


of DERIV4. 

FOR 

using FTLT.FOR with one 

of 

the 

three 

filter 

functions 

after choosing its 

size; 


step 

four- 

plus-i 

transform 

this 

result to obtain 

the 

X, 

second xy, 

or y 

derivative 

of 

the original data 

using 

DERIV4. 

FOR. 





My**- 


The data files used by DERI V4. FOR, 0N0I3E.F0R, 
PILT.FOR, are all unformatted binary random access files. 
The size of the two dimensional input and output arrays 
necessitated the use of such files. Disc I/O time was 
reduced drastically and thus execution time was greatly cut. 
The record length of these data flies is equal to the number 
of data points. 

DERIV4«F0R is the FORTRAN program that performs the 
transformation of the x-derivative of the data, second 
xy-derivative of the data, y-derivative of the data, or the 
transform of the data. The run-tirae parameters ares 

1 • The size of the matrices that hold the input and output 
data- the limit is 64 (for 64 X 64 matrices). Must input an 
integer that is an integral power of two. 

2. The operation to be performed. Enter 1 for the 
transform of x-derivative of the input, 2 for the transform 
of the second xy-derivative of the input, 3 for the 
transform of the y-derivative of the input, 4 for the 
transform of the input. 

3. The type of data- real or complex. Enter 1 for real 
data, a 2 for complex data. 


4. The sign of the transform- minus-i or plus-i. 


Enter -1 


for the ralnus-i transform, 1 for tho olus-i transform. 


5. The output file unit numbers. Enter the unit number for 
the real part data file first, then the unit number for the 
Imaginary part data file. 


6. The Input file unit number(n). Enter the unit number for 
the real part. If the data is complex follow with the unit 
number for the Imaginary part. 
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DERI V4. FOR 

DIMENSION A1 ( 64 , 64 ) , B1 ( 64 , 64 ) , DATA1 (128) 

DIMENSION 01(64,64) 

TYPE 556 

556 FORMAT ( * ENTER SIZE OP MATRIX, LIMIT IS 64 • ,$) 
ACCEPT 557, IZ1 

557 FORMAT ( I ) 

TYPE 558 

558 FORMAT ( ' ENTER D/DU, D2/DUDV , D/DV, OR F(X) ' ,$) 
ACCEPT 557, IZ2 

TYPE 559 

51 .. FORMAT ( ' ENTER REAL( 1 ) OR C0MPLEX(2) DATA TYPE ',$) 
ACCEPT 557, IZ3 
TYPE 560 

560 FORMAT*' ENTER SION OF TRANSFORM (1, -1) ' ,$) 

ACCEPT *>57, IZ5 
IZ4«2*IZ1 

563 FORMAT (21) 

CALL DERIV(IZ1 , IZ2, IZ3, IB4, IZ5, A1 ,B1 , DATA1 ,C1 ) 

STOP 

END 

SUBROUTINE DERIV( INI , IK, IL1 , IN3, ISN1 , A, B, DATA, Cl ) 
DIMENSION A ( IN 1 , INI ) , DATA( IN3 ) ,B( INI , INI ) 

DIMENSION Cl (INI , INI ) 

LOGICAL FLA0,FLA02 
48 FORMAT (21) 

IN2.IN1/2 


IL2-IN1 *IN1 
TYPE 562 

562 FORMAT ( * ENTER OUTPUT FILE ft S(RE,IM) • ,$) 

ACCEPT 48,IF1,IF2 
FLAG2* .TRUE. 

FLAG.. TRUE. 

14 FORMAT (I) 

GO TO (220,230) IL1 
220 TYPE 565 

565 FORMAT (• ENTER FILE H ' ,*) 

ACCEPT 14,IFC 

CALL DEFINE FILE(IFC, IL2,L0C3, 0,0,0) 

666 READ (IFC01) A 

667 CALL REFRMT ( IN 1 , IN2 , A , C 1 ) 

GO TO 333 

230 TYPE 566 

566 FORMAT ( * ENTER FILE ft S(RE,IM)',$) 

ACCEPT 48 i IFA, IFB 

CALL DEFINE FILE ( IFA, IL2, LOCI ,0,0,0) 

CALL DEFINE PTLE ( IFB, IL2, L0C2 ,0,0, 0) 

READ ( IFA^I ) A 
READ (IFB/SM) B 
CALL REFRMT (INI ,IN2,A,C1 ) 

CALL REFRMT ( INI , IN2,B,C1 ) 

333 DO 1 1*1 , INI 

DO 3 II *1 ,IN3 
3 DATA( II )®0 


II IP (FLAG) GOTO 100 
DO 110 K*1 , INI , 2 
A(I,K)— A(I,K) 

110 CONTINUE 

GOTO 111 

100 DO 101 K«2, INI ,2 

A ( I ,K)»-A( I , K) 

mi CONTINUE 

III FLAG*. NOT. FLAG 
DO 4 J=1 , INI 

4 DATA(2*J-1 )=A(I,J) 

IF ( IL1 . EQ. 1 ) GO TO 240 
IF ( FLAG2 ) GOTO 200 
DO 210 K*1 , INI ,2 
B( I , K)*-B( I , K) 

210 CONTINUE 
GOTO 211 

200 DO 201 K=2, INI , 2 
B( I ,K)=-B( I , K) 

201 CONTINUE 

211 FLAG2= . NOT . FLAG2 
DO 241 Jssl , INI 

241 DATA( 2*J )=B( I , J) 

240 CALL MRKFFT( INI ,I3N1 ,IN3, DATA) 

DO 5 J=1 ,IN1 
A(I,J)=DATA(2*J-1 ) 

5 B( I , J)=DATA(2 # J ) 


77 


1 CONTINUE 

DO 60 12*1 , INI-1 
DO 60 J*I2+1 , INI 
TMP*A( 12, J) 

A(I2,J)*A(J,I2) 

A(J,I2)*TMP 

60 CONTINUE 

DO 61 13*1, INI -1 
DO 61 J*I3+1 , INI 
TMP*B( 13, J) 

B(I3.J)*B(J,I3) 

B(J,I3)=TMP 

61 CONTINUE 

DO 2 1=1, INI 
DO 6 J=1 ,IN1 
DATA(2*J-1 )=A( I , J) 

6 DATA(2*J)=B(I,J) 

CALL MRKPPT( INI ,ISN1 ,IN3, DATA) 

DO 7 J=1 ,IN1 

A( I , J)=DATA(2*J-1 ) 

7 B( I , J )=DATA(2*J ) 

1 2 CONTINUE 

2 CONTINUE 

DO 70 12=1, INI -1 
DO 70 J=I2+1 , INI 
TMP=A( 12, J) 

A( 12, J)=A( J , 12) 
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A(J,I2)*TMP 

70 CONTINUE 

DO 71 13*1, INI -1 
DO 71 J*I3+1 , INI 
TMP*B(I3,J) 

B(I3,J)*B(J,I3) 

B(J,I3)*TMP 

71 CONTINUE 

00 T0( 49, 39 ,59, 777) IK 
49 DO 41 1*2, IN2f1 

DO 41 J=IN2+1 , INI 

A(I t J)«A( I ,J)*6 .28318531 # ( J-(IN2+1 ) ) 

41 B(I,J)»B(I,J)*6.28318531*(«J-(IN2+1 )) 
DO 42 I»2, IN2+1 

DO 42 J=2, IN2 

A( I, J)»A( I ,J)*6. 28318531 *(J-( IN2+1 ) ) 

42 B( I, J)=B( I , J)*6. 2831 8531 *( J-( IN2+1 )) 
DO 43 I»IN2+2 , INI 

DO 43 J*2, IN2+1 

A(I,J)sA(I,J)*6.28313531*(J-(IN2+1 )) 

43 B(I, J)=B(I, J)*6.28318531 *(J-(IN2+1 )) 
DO 44 I=IN2+2 , INI 

DO 44 J=IN2+2 , INI 

A( I, J)=A( I , J)*6.2831 8531 *( J-( IN2+1 ) ) 

44 B( I , J)*B( I , J)*6 . 2831 8531 *( J-( IN2+1 )) 
A(1 ,1 )=A(1 ,1 )*(-(IN2+1 ))*1 .57079633 
B(1 ,1 )-B(1 ,1 )*(-(IN2+1 ))*1 .57079633 


mmmmrnwmPMKm 


DO 68 J*2 , INI 

A(1 ,J)*A(l ,J)*(J-(IN2>1 ))*3. 14159265 
B(1 ,J)*B(1 , J)*( J-( IN2+1 ))*3. 14159265 
DO 69 1*2, INI 

A( I , 1 )*A(I ,1 )*3.14159265*(-(IN2+1 )) 

B( 1 , 1 )*B( 1 , 1 )*3.14159265*(-(IN2+1 )) 

GO TO 777 
DO 31 1*2, IN2+1 
DO 31 J*IN2+1 , INI 
A(I,J)*A(I,J)*6.28318531 * 

(J-( IN2+1 ) )*( ABS( I-( IN2+1 ) ) ) 

B( I , J)=B( I , J) # 6 . 28318531 *( ABS( I-( IN2+1 ))) 

DO 32 1*2, IN2+1 

DO 32 J=2,IN2 

A( I , J)*A( I ,J)*6 .2831 8531 * 

(J-(IN2+1 ) )*( ABS( I-( IN2+1 ) ) ) 
B(I,J)*B(I,J)»6.28318531 
*(J-(IN2+1 ))*(ABS(I-(IN2-M ))) 

DO 33 I*IN2+2, INI 
DO 33 J=2, IN2+1 

A(l, J)=A(l ,J)*6 .28^1 8531 *( J-(IN2+1 ))*((IN2+1 )-l) 
B(l,J)=B(l,J)*6. 28318531 *(J-(IN2-M ))*( (IN2+1 )-l) 
DO 34 I=IN2+2, INI 
DO 34 J=IN2+2» INI 

A(I,J)=A(I,J)*6.28318531*('J-(IN2+1 ))*((IN2+1 )-l) 
B(I,J)=B(I,J)*6.28318531*(J-(IN2+1 ))*((IN2+1 )-l) 


- * * -w. 


r^FT*'' 



• 'pm* 


; **. A r+ 


is. ****** T* **" 
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A( 1 , J)=A( 1 , J)*3. 141 59265* ( IN2+1 )»(J-(IN2+1 )) 
67 ii( 1 , J)*B( 1 , J)*3 . 1 41 59265* ( IN2+1 )*(J-(IN2+1 )) 

DO 66 1*2, INI 

A(I,1 )*A( I , 1 )*3 • 1 41 59265*( I-( IN2+1 ))*(IN2+1 ) 
66 B(I,1 )*B(I,1 )*3 • 1 41 59265* ( I-( IN2+1 ) )*( IN2+1 ) 

A(1 ,1 )*A(1 ,1 )*(IN2+1 )*( -( IN2+1 ))*1 .57079633 
B(1 ,1 )*B( 1 ,1 )*( IN2+1 )*( -( IN2+1 ))*1 .57079633 
00 TO 777 

59 DO 51 1*2 , IN2+1 

DO 51 J*IN2+1,IN1 

A( I , J)*A( I ,J)*6 .2831 8531 *( ABS( I-( IN2+1 ))) 

51 B(I,J)=B(I,J)*6. 2831 8531 *( ABS( I-( IN2+1 ))) 

DO 52 1*2 , IN2+1 

DO 52 J=2 , IN2 

A( X , J)=A( I , J)*6 . 28318531 *ABS( I-( IN2+1 )) 

52 B( I , J)=B( I , J)*6. 2331 8531 *( ABS( I-( IN2+1 ))) 

DO 53 I=IN2+2 , TNI 

DO 53 J=2 , IN2+1 

A(I,J)=A(I, J)*6. 2831 8531 *( (TN2+1 )-l) 

53 B(I,J)=B(I, J)*6. 2831 8531 *((IN2+1 )-I) 

DO 54 I=IN2+2 , INI 

DO 54 J=IN2+2, INI 

A(I,J)=A(I,J)*6. 28318531 *( (IN2+1 )-I) 

54 B(I, J)=B(I,J)*6.2831S531*( ( IN2+1 )-I) 

A( 1 ,1 )=A(1 ,1 )*(IN2+1 )*1 .57079633 

B(1 ,1 )=B(1 ,1 )*( IN2+1 )*1 .57079633 
DO 65 J=2 , INI 


A(1 ,J)«A(1 ,J)*3.14159265*((IN2+1 )) 

65 B(l ,J)*3.14159265*(IN2+1 ) 

DO 64 1*2, INI 

A( 1 , 1 )*A( 1 , 1 )*3.14159625*(I-(IN24l )) 

64 B(I,1 )*B(I,1 )*3.14159625 # (I-(IN2+1 )) 

C THE REAL PART OP THE TRANSFORM IS IN 21 

C THE IMAO. PART OP THE TRANSFORM IS IN 22 

777 CALL DEFINE PILE( IP1 , IL2,L0C4, 0,0,0) 

CALL DEFINE PILE(IP2, IL2,LOC5, 0,0,0) 
WRITE (IF201) A 
WRIT'5 ( IP1 #1 ) B 
700 HETUrtN 

END 

SUBROUTINE MPKFFT(NN,ISIGN,IQ1 , DATA) 
DIMENSION DATA(IQI) 

N« 2*NN 
C 
C 
C 
C 

C PAST FOURIER TRANSFORM ROUTINE 

C 

C 

J*1 

DO 5 1*1 ,N,2 
IF( I-J ) 1 ,2,2 
1 TEMPR=DATA( J) 


TEMPI«DATA( J+1 ) 

DATA(J)-DATA(I) 

DATA( J+1 )*DATA( 1+1 ) 

DATA( I )«TEMPR 
DATA( 1+1 )*TEMPI 
M.N/2 

IF( J-M)5,5,4 

J*J-M 

M*M/2 

IF(M-2)5,3,3 

J*J+M 

MMAX.2 

IF(MMAX-*N)7, 10, 10 
ISTEP*2*MMAX 

THETA«6. 2831 853 /FLOAT ( I SIGM*MMAX) 
SINTH=8IN ( THETA/2 ) 
WSTPR»-2*SINTH*SINTH 
WSTPI sSIN ( THETA ) 

WR=1 

WI=0 

DO 9 M*1 ,MMAX, 2 
DO 8 I=M, N, ISTEP 
Jal+MMAX 

TEMPR=WR*DATA(J)-WI*DATA(J+1 ) 
TEMPI*WR # DATA( J+1 )+WI*DATA(<J) 

DATA ( J ) =DAT A ( I ) -TEM PR 
DATA( J+1 )=DATA( 1+1 ) -TEMPI 


DATA ( I ) *D AT A ( I ) +TEMPR 
DATA(I+1 )«DATA(I+1 )+TEMPl 
TBMPR.WR 

WR.WR* W3T PR-W I * WST P I +WR 
WI.WI*WSTPR+TEMPR*WSTPI+WI 
MMAX.ISTBP 
00 TO 6 


RETURN 


SUBROUTINE REFRMT(IS,IS2,A,C) 
DIMENSION C( IS2, IS2) t A( IS, IS) 
FORMAT ( I ) 

DO 60 1*1 ,IS2 
DO 60 J*1 , IS2 
C( I , J)*A( I , J) 

DO 61 I»IS2+1 , IS 
DO 61 J=IS2+1 , IS 
A(I-IS2,J-IS2)*A(I,J) 

A(I , J)*C( I-IS2, J-IS2) 

DO 62 1*1 ,IS2 
DO 62 J=1 ,IS2 
C(I,J)=A( I , J+IS2) 

DO 63 1=1 ,IS2 
DO 63 J=1 , IS2 
A(I,J+IS2)=A(I+IS2,J) 

A( I+IS2, J)=C( I , J) 





FORMAT (8F) 
RETURN 
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PILT.POR is the PORTRAN filter program. The choice of 
filter functions is ideal square and circular low-pass 
filter functions, and a flat- topped filter function. The 
extent of the filter function can be varied by changing the 
value of M for the square and flat-topped filter function 
and R (radius) for the circular filter function. The 
run-time parameters are: 

1. The size of the square matrices- same conditions hold as 
in DERI V4. FOR. 

2. The input and output files. The data type is assumed to 
be complex, therefore a total of four logical unit numbers 
must be entered. The input unit numbers must be first, and 
for each pair of unit numbers the real unit number is first. 

3. The filter function desired. Enter 1 for the square 
filter function, 2 for the circular filter function, 3 for 
the flat-topped pyramid filter function. 

4. The size or extent of the filter function- limit is 
one-half the size of the matrices used by the program to 


hold data. 
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&$-$?- }**-**» m : 


FILT.POR 

DIMENSION A1(64,64),B1(64,64) 

TYPE 10 

10 PORMATC ENTER SIZE OP MATRIX ' ,$) 

ACCEPT 11'ISZI 

11 PORMAT ( I ) 

ISZ2.ISZ1/2+1 

CALL PILTR(A1 ,B1 , ISZ1 , ISZ2) 

STOP 

END 

SUBROUTINE PILTR( A,B, INI , IORIO) 

DIMENSION A( INI , INI ) ,B( INI , INI ) 

IL2vIN1 # IN1 

401 PORMAT ( I ) 

402 PORMAT (21) 

TYPE 460 

460 PORMAT ( * ENTER INPILES , OUTFILES ( RE , IM ) •,$) 

ACCEPT 461 ,IF1 ,IF2,IF3,IF4 

461 PORMAT (41) 

CALL DEFINE PILE ( IP1 ,IL2, LOCI ,0,0,0) 

CALL DEFINE FILE(IF2,IL2,L0C2, 0,0,0) 

CALL DEFINE FILE(IF3, IL2,L0C3, 0,0,0) 

CALL DEFINE FILE(IF4, IL2,L0C4,0, 0,0) 

READ(lF1/h) A 
READ(IF2#1 ) B 
TYPE 410 

410 PORMATC IS FILTER SQ( 1 ) ,CIRC(2 ) ,CUT PYRAM(3) ’ ,$) 



« 


I 

1 




i 



\ 


I 


\ 


Mix 


ACCEPT 401 , ISHAPE 

00 TO (420,430,450) ISHAPE 

420 TYPE 421 

421 PORMAT ( * ENTER M( SI ZB-2*M+1 ) \$) 
ACCEPT 401 , IMM2 

00 TO 440 

430 TYPE 431 

431 PORMAT (' ENTER RAD ' ,$) 

ACCEPT 401.IRAD 

440 IP ( ISHAPE. EQ i 1 ) IMID-IMM2 
IP (ISHAPE. BQ. 2) IMID-IRAD 
DO 441 I«1 ,I0RI0-IMID-1 

DO 441 J«1 ,IN1 
A(I, J)«0.0 
B( I , J)»0.0 

441 CONTINUE 

DO 442 I.IORIG+IMID+1 ,IN1 
DO 442 J-1 ,IN1 
A(I,J)-0.0 
B(X,J)-0.0 

442 CONTINUE 

DO 443 I*I0RI0-IMID, IORIG+IMID 
DO 443 J*1 ,I0RI0-IMID-1 
A(I,J)-0.0 
B(I,J)*0.0 

443 CONTINUE 

DO 444 I-IORIG-IMID, IORIG+IMID 


DO 444 JaIORXQ+IMID+1 , INI 

A(l,J)-0.0 

B(I # J)«0.0 

CONTINUE 

IF (I8HAPE.EQ.1 ) 00 TO 490 

K1 »IN1 /2+TMID+2 

J3-(INl/2-IMID)*2 

K2.IN1/2-IMID+1 

J4.J3+4*IMID+2 

K3-K1 -1 

J5.J3+1 

J6-J4/2 

RIGNalORIO 

DO 140 I2.I0RI0-IMID,I0RI0+IMID 

DO 140 I3»I0RIG-IMID , IORIO+IMID 

L1-I2-IN1/2-1 

L2«(I3+1 )/2-IN1/2-1 

XR«(I2-RI0N)**2 

YR«(I3-RI0N)**2 

R-IMID 

IF (SQRT(XR+YR).LE.R) GOTO HO 

A(I2,I3)»0.0 

B(I2, I3)»0.0 

CONTINUE 

00 TO 490 

CIRCLE FILTER IS FINISHED 


450 TYPE 451 

451 FORMAT ( ' ENTER MID(PLATBAU«2*MID+1 ) ' ,$) 
ACCEPT 401 , IMID 

IRMP»( IORIG-1 )-( IMID+I ) 

XDBL.1 .0 

DO 22 I*3,I0RI0-IMID-1 
DO 23 J«I,IN1-(I-2) 
A(I,J)-(XDBL/IRMP)*A(I,J) 
B(I,J)-(XDEL/IRMP)*B(I,J) 

23 CONTINUE 
XDEL«XDEL+1 .0 

22 CONTINUE 

XDEL-1 .0 
K-0 

DO 24 I-IN1-1 .I0RI0+IMID+1 ,-1 
DO 25 J»3+K, INI-1 -K 
A(I,J)*(XDEL/IRMP) # A(I,J) 
B(I,J)«(XDEL/IRMP)*B(I,J) 

25 CONTINUE 

XDEL-XDEL+1 .0 
K.K+1 

24 CONTINUE 
XDEL-1 .0 

DO 26 J*3# IORIG-IMID-1 
DO 27 I*J+1 , IN1-J+1 

A( I , J)» (xd.bl/irmp)*a( I , j) 

B(l, J)*(XDEL/IRMP) # B(I , J) 


27 CONTINUE 
XDBL-XDBL+1 .0 

26 CONTINUE 

XDEL-1 .0 
K*0 

DO 28 J-IN1-1 , IORIO+IMID+1 
DO 29 I*4+K, IN1-2-K 
A( I , J)*(XDEjj/lRMP)*A( I , J) 
B( I, J)*(XDEL/lRMP)*B( I, J) 

29 CONTINUE 
K-K+1 

XDEL«XDEL+1 .0 

28 CONTINUE 

DO 30 J*1 , INI 
A(IN1 ,J)=0.0 
B( INI ,J)«0.0 

30 CONTINUE 

DO 31 1*1, INI 
A( I , INI )»0.0 
B(I ,IN1 )*0.0 

31 CONTINUE 

DO 20 1*1,2 
DO 20 J=1 , INI 
A(I,J)=0.0 
B(l , J)*0.0 
20 CONTINUE 




adds Gaussian 
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QNOISE.FOR is the FORTRAN program that adds Gaussian 
ord inant dependent noise to the input data. The input data 
is not destroyed. The run-time parameters are: 


1 . The sise of the matrices- the same conditions hold as for 
DERI V4. FOR. 


2. The scale factor of the noise. 


?. The type of data- real or complex. Enter 1 for real, 2 
for complex. 

4. The input logical unit number(s), the unit number(s) that 
is (are) to contain signal plus noise, the unit number(s) 
that is (are) to contain the noise only. If the data is 
complex two unit numbers are required for each case, with 
the real unit number being the first in all cases. 



f 


w 
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GN0I3B.P0R 

DIMENSION A(64.64),B(64,64) 

REAL S? 

TYPE 10 

10 PORMAT ( ' ENTER SIZE OP MATRIX ’ ,$) 

ACCEPT 11,ISZ 

11 PORMAT ( I ) 

TYPE 12 

12 FORMAT (• ENTER SCALE FACTOR • ,$) 

ACCEPT 1 3,SF 

13 PORMAT (0) 

CALL NOIS( A,B,SF, ISZ) 

STOP 

END 

SUBROUTINE NOIS( A1 ,B1 , SF1 , ISZ1 ) 

DIMENSION A1 ( ISZ1 , ISZ1 ) ,B1 ( ISZ1 ,ISZ1 ) 

REAL 3P1 
IL2=ISZ1 *ISZ1 
K0UNT*0 
TYPE 100 

100 PORMATC REAL(1) OR COM PLEX(2) DATA TYPE ' ,$) 

ACCEPT 110, IK1 

111 PORMAT (21) 

110 PORMAT ( I ) 

112 PORMAT (2G) 

PORMAT (G) 

PORMAT (31) 


Jj 

i 

*! 


a 



113 

114 
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115 FORMAT (61) 

00 T0(1 20, 130) IK1 

120 TYPE 121 

121 FORMAT (' ENTER INFILE, S+N,N PILE • ,*) 

ACCEPT 1 14, IF1 , IF81 , IFNI 

CALL DEFINE FILE(IF1 ,IL2, LOCI ,0,0,0) 
CALL DEFINE FILE( IFS1 , IL2, L0C2, 0,0,0) 
CALL DEFINE FILE( IFNI , IL2, L0C3, 0,0,0) 
READ(IF1/?1 ) A1 
00 TO 200 

130 TYPE 131 

131 FORMAT (' ENTER INPUT RE,IM,(S+N) 

! RE, IM, (N) RE, IM' ,$) 

ACCEPT 115* IP1 ,IF2,IFS1 ,IFS2,IFN1 ,IFN2 
CALL DEFINE FILE(IF1 ,IL2, LOCI ,0,0,0) 
CALL DEFINE FILE(IF2, IL2,L0C2, 0,0,0) 
CALL DEFINE FILE(IFS1 , IL2,LOC3, 0,0,0) 
CALL DEFINE FILE(IFS2, IL2, L0C4, 0,0,0) 
CALL DEFINE FILE(IFN1 , IL2, L0C5, 0,0,0) 
CALL DEFINE FILE( IFN2, IL2, L0C6, 0,0,0) 
READ (IF1#1) A1 

200 RMS *0.0 

KOUNT.KOUNT+1 
DO 300 I*1,IS551 

DO 300 J=1 , ISE1 

K*. 

P*RAN ( 5 ) 


S.RAN(IO) 

IP (A1 (I, J).LT.1 .0E-20) XN*2.*P # A1 (I, J) 
IP (A1 (I,J).LT.1 .0E-20) 00 TO 250 
XN»ABS(A1(I,J))* 

! (-AL00(P*SQRT(6.28318*SP1*ABS(A1 (I,J))))) 
XN*SQRT(2.0*3F1*XN) 

IP (S.0T.0.5) XN— XN 
TEMP.A1 (I,J) 

A1(I,J)=A1(I , J)+XN 
B1 (I,J)*XN 

IP (A1 (I.J).LT.O.O) A1 ( I , J)*0.0 
RMS*(A1 (I , J)-TEMP) # *2+RMS 
CONTINUE 

RMS*SQRT(RMS/IL2) 

TYPE 1 1 3 > RM8 

IP (KOUNT .EQ. 2) GOTO 500 
WRITE (iFSI/jM) A1 
WRITE (IFNI/Pl) B1 
IP ( IK1 .EQ. 1) GOTO 400 
IP (KOUNT .EQ. 1) READ (lF2/h) A1 
IP (KOUNT .EQ. 1 ) GOTO 200 
WRITE (IFS2/SM) A1 
WRITE (IFN2/5M) B1 


RETURN 
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